use ${newdata}carownership_dataset_bergen, clear

drop if $trmgroup == . | $trmgroup == 0
capt drop bergen 
capt drop aux 

egen long famid_num = group(familienr)

gen commuter = ($trmgroup == 1 | $trmgroup == 3)
gen bergen = ($trmgroup == 1 | $trmgroup == 2)

local year_vars
local beta_int
local year_grk
forvalues i = $firstyearpre/$lastyearpost {
	if `i' != $lastyearpre qui gen year_`i' = (year == `i')
	if `i' != $lastyearpre local beta_int `beta_int' c.year_`i'#c.commuter
	if `i' != $lastyearpre local year_vars `year_vars' c.year_`i'
	if `i' != $lastyearpre local year_grk `year_grk' c.year_`i'#i.grk	
}	

/* All time-varying variables based on geography (residence or work location)
   are set to their 2014 values (variation post 2014 is endogenous to treatment) */
foreach var in dist time_work PublicVSCarTime_fam_mean PublicDiffCarTime_fam_mean grk {
	replace `var' = . if year != 2014
	bysort familienr (`var'): replace `var' = `var'[_n-1] if missing(`var')
}	

/* === REGRESSION FOR BERGEN =================================================*/
reghdfe $yvar /// LHS
	`beta_int' /// Years*treatment (base == `baseyear')
	`yearvars' /// Years (base == `baseyear')
	$xvar /// other controls
	commuter /// Absorbs the difference in the base year
	if ($trmgroup == 1 | $trmgroup == 2) /// paying + non-paying in Bergen
		& year >= $firstyearpre /// 
		& year <= $lastyearpost /// 
	, absorb( ///
		i.grk /// Grunnkrets FE
		`year_grk' /// Year*grunnkrets FE (base == `baseyear')
		i.famid_num /// Family FE
	) ///
	vce(cluster i.grk)
estimates save ${ster}figC1_${yvar}_1, replace
/* Saving estimates */
matrix b = e(b)
matrix v = e(V)

/* Storing them as variables */
qui gen beta = .
qui gen se = .
qui gen xaxis = .
local coef = $lastyearpost - $firstyearpre /* The number of coef we want to keep */
forvalues i = 1/`coef' {
	qui replace xaxis = $firstyearpre + `i' - 1 in `i'
	qui replace beta  = b[1,`i'] in `i' // the i-th coefficient of b
	qui replace se    = sqrt(v[`i',`i']) in `i' // the i-th diagonal element of v
}
/* Adding the base year here */
qui replace xaxis = xaxis + 1 if xaxis >= $lastyearpre /* To skip the base year */	
local lastobs = `coef' + 1 /* Adding an observation for the base year */
qui replace xaxis = $lastyearpre 	in `lastobs'
qui replace beta = 	0 				in `lastobs' /* Coef is normalized to zero */
qui replace se = 	0 				in `lastobs'
/* 95 % CI */
qui gen plus2se = beta + 1.96 * se
qui gen minus2se = beta - 1.96 * se

/* === REGRESSION FOR OTHER CITIES ===========================================*/

reghdfe $yvar /// LHS
	`beta_int' /// Years*treatment (base == `baseyear')
	`yearvars' /// Years (base == `baseyear')
	$xvar /// other controls
	commuter /// Absorbs the difference in the base year
	if ($trmgroup == 3 | $trmgroup == 4) /// paying + non-paying in other cities
		& year >= $firstyearpre /// 
		& year <= $lastyearpost /// 
	, absorb( ///
		i.grk /// Grunnkrets FE
		`year_grk' /// Year*grunnkrets FE (base == `baseyear')
		i.famid_num /// Family FE
	) ///
	vce(cluster i.grk)
estimates save ${ster}figC1_${yvar}_2, replace

estimates use ${ster}figC1_${yvar}_2
/* Saving estimates */
matrix b2 = e(b)
matrix v2 = e(V)

/* Storing them as variables */
qui gen beta2 = .
qui gen se2 = .
qui gen xaxis2 = .
local coef = $lastyearpost - $firstyearpre /* The number of coef we want to keep */
forvalues i = 1/`coef' {
	qui replace xaxis2 = $firstyearpre + `i' - 1 in `i'
	qui replace beta2  = b2[1,`i'] in `i' // the i-th coefficient of b
	qui replace se2    = sqrt(v2[`i',`i']) in `i' // the i-th diagonal element of v
}
/* Adding the base year here */
qui replace xaxis2 = xaxis2 + 1 if xaxis2 >= $lastyearpre /* To skip the base year */	
local lastobs = `coef' + 1 /* Adding an observation for the base year */
qui replace xaxis2 = $lastyearpre 	in `lastobs'
qui replace beta2 = 	0 			in `lastobs' /* Coef is normalized to zero */
qui replace se2 = 	0 				in `lastobs'
/* 95 % CI */
qui gen plus2se2  = beta2 + 1.96 * se2
qui gen minus2se2 = beta2 - 1.96 * se2

/* === FIGURE ============================================================*/

/* Time dimension */
gen time = mdy(1,1,xaxis )
format time %td

gen ttime = year(time)

local announcement_date   = 2014 + (mdy(2,18,2015) - mdy(1,1,2015) + 1) / 365
local implementation_date = 2015 + (mdy(2,1,2016) - mdy(1,1,2016) + 1) / 366

twoway ///
	(rarea plus2se2 minus2se2 ttime ///
		, sort lcolor(gs10%50) fcolor(gs10%50)) ///	
	(rarea plus2se minus2se ttime ///
		, sort lcolor(gs5%50) fcolor(gs5%50)) ///
	(scatter beta ttime ///
		, sort connect(l) lcolor(black) mcolor(black) lpattern(solid)) ///
	(scatter beta2 ttime ///
		, sort connect(l) lcolor(black) mcolor(black) lpattern(dash)) ///	
	, graphregion(color(white))	///
	xtitle("") ytitle("Diff. betw. paying and non-paying commuter") ///
	name(pretrend_triple, replace) ///
	legend(on order(3 4) ///
		label(3 "Estimate, Bergen") ///
		label(4 "Estimate, control cities") ///
		pos(6) col(2) ///
	) ///
	xline(`announcement_date', lpattern(shortdash) lcolor(black)) /// mdy(2,1,2016) - implementation date
	xline(`implementation_date', lpattern(dash) lcolor(black)) /// mdy(2,18,2015) - announcement date
	yline(0) xlabel(2011(1)2017) ylabel(-0.08(0.04)0.08)
	
	
graph export "${figures}figC1${yvar}.png", as(png) replace
graph export "${figures}figC1${yvar}.pdf", as(pdf) replace
